 	PROGRAM COLLISIONMATRIX
 
C********************************************************************************
C	NB I have included an analytic expression for Ar elastic under 1 eV.
C	Check to see if appropriate
C********************************************************************************

C     Code calculates collision matrix for binary gas mixtures
C	gas 1 undergoes elastic & inelastic collisions & ionization & attachment 
C	gas 2 undergoes elastic & inelastic collisions & ionization & attachment

c	OUTPUT: 
c       (i)  Dimensionless a's in the arrays A(nu,nud,l)
c	  (ii) Dimensionlesss interaction integrals (U(i,j)'s) in vrate(nq,nu)
c		   required for the rate calculations.
c	N.B. Collision matrix and interaction integrals for mixtures are
c	made dimensionless with respect to gas 1. This is tested for in DCCONS.

c	Files output from COLLISIONMATRIX to read into DCCONS:
c		(i)  fort.14 = dimensionless collision matrix
c		(ii) fort.13 = data for consistency checks
c	    (iii)fort.12 = quantities for calculation of rate coefficients

c	INPUT:

c	(i) Cross-sections are read in as arrays where J identifies 
c         the particular processes.

c	NO(J) = no. of data points in the cross-sections
c	NP(J) = no. of partial cross-sections for jth process
c	EN(J,K) = kth energy point for the jth cross-section
c	QS(J,K,L) = Lth partial-cross-section evaluated at en(j,k)
c	SW(J) = n_o(lower)/n_o = fractional population in the lower state of the
c		  jth process
c       EN(J,2) = threshold for the jth collisional process
c	NQ1 = no. of cross-sections of gas 1
c	NQ2 = no. of cross-sections of gas 2
c	NQI1 = no. of non-elastic processes in gas 1
c	NQI2 = no. of non-elastic processes in gas 2
c	CON1 = conc. of gas 1

c	(ii) Input data

c	TI = Basis temperature
c	T = Neutral temperature
c	W1 = mass of neutral gas 1 element (NB. Boltz solver must be scaled to this value)
c	W2 = mass of neutral gas 2 element
c	W = mass of the charged particle (generally set to electron mass)

c	(iii) Numerical input

c	RMAX = The maximum value of \nu_{max}-1 in inversion routine
c	LMAX = The maximum value of l_{max}-1 in inversion routine
c	NEL = No. of intervals for elastic interaction integrals
c	WEL = The maximum energy range is given by 1/(1-WEL) 
c	NINEL = No. of intervals for non-elastic interaction integrals
c	WINEL = The maximum energy range is given by 1/(1-WEL) 
c	NGQ = No. of guass quadrature points.

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

      CHARACTER datetime*24,lable*80
      INTEGER R,S,RMAX
      COMMON/WEALTH/FACLOG(150),GAMLOG(150),VX(10,100,100),A(100,100,9),
     +   VXLIQ(10,100,100),VXATT(10,100,100),VXELAS(10,100,100),
     +      VXINEL(10,100,100)
      COMMON/COL/VXION(10,100,100) 
      COMMON/GOOD/EN(21,1100),QS(21,1100,0:10),NO(21),npq(21)
      COMMON/LAW/SW(21),DE(21),V(21),VI(21)
      COMMON/INLAW/nq,nq1,nq2,nqi1,nqi2
      COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)

c ---	Data input  ----------------------------------------------------------

	open(unit=3,file='Tb.dat',status='old')
	
!      open(unit=1,file='positron_argon_gas_McEachran_aniso.dat',      
!     &  status='old')
!	 open(unit=1,file='NeHayashi.dat',status='old')
	 open(unit=1,file='Hgmodel.dat',status='old')

!      open(unit=1,file='Water_structured.dat',status='old')
     
      !open(unit=1,file='positron_argon_liquid.dat',status='old')
       READ(1,*)lable		! Label for collision matrix
       READ(1,*)nq1		! No. of xsections for gas 1
       READ(1,*)nq2		! No. of xsections for gas 2
       READ(1,*)W1		! Mass (amu) for gas 1 
       READ(1,*)W2		! Mass (amu) for gas 2
       READ(1,*)LMAX		! Highest order lmax required
       READ(1,*)RMAX		! Highest order numax required
       READ(3,*)TI		! Basis temperature (Tb)
       READ(1,*)T		! Neutral gas temperature
       READ(1,*)NEL,WEL,NINEL,WINEL ! No. of intervals for integrations
		    ! Width of entire integration domain = 1/1-W
	 NGQ = 200
       READ(1,*)CON1		! Concentration of gas 1
       CON2=1.D0-CON1		! Concentration of gas 2
       nq = nq1 + nq2		! Total number of qs
       nqi1 = nq1 - 1		! No. of non-elastic cross-sections gas 1
       nqi2 = nq2 - 1		! No. of non-elastic cross-sections gas 2

       WRITE(2,1005)lable,datetime
       WRITE(2,2504)NQ,T,W1,W2
       WRITE(2,2505)TI,LMAX,RMAX,CON2
	
       TK=8.617385D-5*T				! kT_o
       TKI=8.617385D-5*TI			! kT_b

C      READ TABULATED CROSS SECTIONS

       DO 100 I=1,NQ
       READ(1,*)NO(I),npq(i),de(i),sw(i),lable
       PRINT*, NO(I),npq(i),de(i),sw(i),lable

		! NO(I) = no. of data points for i
	 	! npq(i) = no. of partial cross-sections
	        ! de(i) = threshold for process i
		! lable = label for process i

      ! NB. de(i) must correspond to en(i,2) in the input data

       J=NO(I)
       READ(1,*)(EN(I,K),K=1,J)		! Energy points (eV)
	print*, 'npq(i)=',npq(i)
       do 101 L=0,npq(i)
         READ(1,*)(QS(I,K,L),K=1,J)     ! Partial cross-sections (ang^2)
         WRITE(6,*)(QS(I,K,L),K=1,J)     ! Partial cross-sections (ang^2)
         write(6,*) L
 101    continue	

	 WRITE(10,*) no(i),de(i),sw(i),lable
       WRITE(10,*)
       WRITE(10,1151)(EN(I,K),K=1,J)
       WRITE(10,*)
       do 102 L=0,npq(i)	
        WRITE(10,1151)(QS(I,K,L),K=1,J)
        write(10,*)
102    continue

       V(I)=0.D0
       VI(I)=0.D0
       IF(T.GT.0.D0)V(I)=DEXP(-DE(I)/TK)	
       IF(TI.GT.0.D0)VI(I)=DEXP(-DE(I)/TKI)

 100   CONTINUE

 1151	FORMAT(9F11.4)

       W=5.48579903D-4		! Electron mass in amu
       WM1=W/(W1+W)  		! Calculating mass ratios
       WM2=W/(W2+W)

       LMAX1=LMAX+1	! This accounts for shift in array bounds 0->1
       IRMAX1=RMAX+1

	 DO 11 L=1,LMAX1		! Initialisation
	 DO 11 I=1,IRMAX1
	 DO 11 J=1,IRMAX1
	  A(I,J,L)=0.D0
	  VXELAS(L,I,J)=0.D0
	  VXINEL(L,I,J)=0.d0
	  VXATT(L,I,J)=0.D0
	  vxion(L,I,J)=0.d0
11	  VX(L,I,J)=0.D0

C	MAKE ARRAY FACLOG (LOG OF FACTORIAL OF INTEGERS)

	FACLOG(1)=0.D0
	DO 919 I=2,RMAX+3
 	XI=DFLOAT(I-1)
919	FACLOG(I)=FACLOG(I-1)+DLOG(XI)

C	MAKE ARRAY GAMLOG (LOG OF GAMMA FUNCTION OF HALF INTEGERS)

	GAMLOG(1)=0.572364942924700176D0
	DO 921 I=2,RMAX+LMAX+4
	XI=DFLOAT(2*I-3)/2.D0
921	GAMLOG(I)=GAMLOG(I-1)+DLOG(XI)


C  Calculate Integrals For Gas 1
	
        if(con2.eq.1.d0)go to 223
	NO_GAS=1

	n_elast1=nq1          ! Process number for elastic process gas 1
        CALL VINT(LMAX1,IRMAX1,TKI,NGQ,NEL,WEL,n_elast1,NO_GAS)
		! VINT calculates the interaction integrals for elastic
		! collision processes, which are stored in the process
		! number n_elastic1

	n_inel1_min=1	     ! Process number of first inelastic for gas 1
	n_inel1_max=nq1-1
       CALL UINT(LMAX1,IRMAX1,TKI,TK,NGQ,NINEL,WINEL,
     & 					n_inel1_min,n_inel1_max,NO_GAS)
		! UINT calculates the interaction integrals for non-elastic
		! conservative collision processes, which are stored in 
	        ! the process numbers n_inel through to nq1-1
		
	 n_ion=nq1-1         ! Process number for elastic process gas 1
 !       CALL RINT(LMAX1,IRMAX1,TKI,NGQ,1999,WINEL,0.5d0,1,n_ion,NO_GAS)
		! RINT calculates the interaction integrals for 
		! ionization collision processes, which are stored in 
	        ! the process number n_ion

!	 n_att=nq1-1          ! Process number for elastic process gas 1
!        CALL ATTINT(LMAX1,IRMAX1,TKI,NGQ,NINEL,WINEL,n_att,NO_GAS)
		! ATTINT calculates the interaction integrals for 
		! attachment collision processes, which are stored in 
	        ! the process number n_att

        CALL COLMAT(WM1,CON1,LMAX,RMAX,T,TI,NO_GAS)
		! COLLMAT converts the dimensionless interaction integrals
		! into the dimensional big V's


C  Calculate Integrals For Gas 2

      IF(CON2.EQ.0.D0)GOTO 222

	
c	Reset VX to zero for second gas calculations.
	
	do 105 l=1,lmax1
	do 105 i=1,irmax1
	do 105 j=1,irmax1
	  vxion(l,i,j)=0.d0
 105	  vx(l,i,j) = 0.d0
	
223	continue	

	NO_GAS=2
      CALL VINT(LMAX1,IRMAX1,TKI,NGQ,NEL,WEL,nq,NO_GAS)
 !     CALL UINT(LMAX1,IRMAX1,TKI,TK,NGQ,NINEL,WINEL,nq1+1,nq-2,NO_GAS)
 !     CALL RINT(LMAX1,IRMAX1,TKI,NGQ,NINEL,WINEL,0.5d0,1,nq-1,NO_GAS)
 !     CALL ATTINT(LMAX1,IRMAX1,TKI,NGQ,NINEL,WINEL,nq-1,NO_GAS)

      CCC=WM1*W1/(WM2*W2)
      CCC=CON2*DSQRT(CCC)		
 		! The collision matrix output is the dimensionless matrix
		! i.e. dsqrt(2kT_b/mu1)*sigma0 is not included
	        ! where mu1=reduced mass of electron and neutral 1
		! This above line makes sure that the total collision matrix for the
		! mixture are scaled with respect to the reduced mass 1
       
	CALL COLMAT(WM2,CCC,LMAX,RMAX,T,TI,NO_GAS)

222   CONTINUE

	DO 789 L=1,LMAX+1 
		! Output the total interaction integral (dimensionless)
	 L1=L-1			! - the ionization interaction integral
	 WRITE(9,49)L1
	 DO 789 I=1,10
	  IK=I
	  IF(L.EQ.1)IK=10
789	  WRITE(9,1234)(VX(L,I,J),J=1,IK)
1234	FORMAT(1X,10D13.5)
49	FORMAT(/3H L=,I2)
	
	! This converts from dimensionless J's to dimensionless a's.
	! i.e. the output is a(nu,nud,l)/[sigma0*dsqrt(2kT_b/mu1)]
	
	DO 2 L=1,LMAX
        DO 2 R=1,RMAX
	 FRL=FACLOG(R)-GAMLOG(R+L)
         DO 2 S=1,RMAX
	  FSRL=FACLOG(S)-GAMLOG(S+L)-FRL
	  FSRL=DEXP(FSRL)
	  A(R,S,L)=DSQRT(FSRL)*A(R,S,L)
2     CONTINUE

	DO 791 L=1,2
	 L1=L-1
	 WRITE(17,49)L1
	 DO 791 I=1,5
 791	  WRITE(17,1234)(A(I,J,L),J=1,5)
	WRITE(14) A
	WRITE(13) LMAX,RMAX,T,TI,W,W1,W2
	WRITE(12) nq1-1,nq2-1,vrate0,vrate1

	
c	Estimate of the E/n_o(Td) for this T_b (based on Maxwell model)

	E=DSQRT(-0.5D0*A(2,1,1)*A(1,1,2))
	EON=E*DSQRT((W1+W)/W1)*TI*1.72341D-3
	WRITE(2,881)EON
	WRITE(6,881)EON
	

c	Format statements for output

 881	FORMAT(' E/N = ',1PD11.4,' Td.')
 1005 FORMAT(a,a)
 1105 FORMAT(/10(1PE12.5))
 2502 FORMAT('0',10X,'ANGULAR DISTRIBUTIONS FOR'/10X,A/)
 2501 FORMAT('0',10X,6(F6.3)/10X,6(F6.3)/10X,6(F6.3))
 2503 FORMAT('0',10X,'ALL INELASTICS MULTIPLIED BY  ',F6.3)
 2504 FORMAT('0',10X,'NUMBER OF CROSS SECTIONS',I4/'0',10X,
     1'GAS TEMPERATURE	 ',F8.3/'0',10X,'MOLECULAR WEIGHTS  ',2F9.4)
 2505 FORMAT('0',10X,'TEMPERATURE PARAMETER	',F9.1/'0',10X,'LMAX =
     1 ',I2,5X,'RMAX = ',I3/'0',10x,'Concentration of gas 2',f15.6)
 2506 FORMAT('0',10X,'STATISTICAL WEIGHTS'/5(10x,6(1pd12.4)/))

c 57	CALL EXIT
	END

c ****************************************************************************

      SUBROUTINE VINT(LLIM,NLIM,TKI,NGQ1,NPOL,WI,JJJ,no_gas)

c 	This subroutine calculates the dimensionless v's for the elastic interactions.
c  	It doesn't calculate the V's or U's.
c	Ensure the jjj cross-section is the elastic cross-section

c	NB. The calculationof vrates for elastic must go from l=0 not l=1
c	as in the traditional rate calculations.

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/WEALTH/FACLOG(150),GAMLOG(150),VX(10,100,100),A(100,100,9),
     +   VXLIQ(10,100,100),VXATT(10,100,100),VXELAS(10,100,100),
     +      VXINEL(10,100,100)
      COMMON/GOOD/EN(21,1100),QS(21,1100,0:10),NO(21),npq(21)
      COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)
      COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)

      
	CALL COTES(NPOL,WI,1)
	SCA2=1.D-26   	! This is the scaling parameter used in cotes

	do n=1,nlim
c	 do l=0,1
	  vrate0(no_gas,0,n-1)=0.d0
	  vrate1(no_gas,0,n-1)=0.d0
c	 enddo
	end do


	DO 11 LL=1,LLIM	! L=0 Interaction integral is zero, but rate
			! calculations require the l=0
	 L=LL-1
	 WRITE(6,505)L
         FL=DFLOAT(L)
	 AL=FL+.5D0
       
       DO 2 N=1,NLIM 	 ! Initialisation
        DO 2 NP=1,NLIM
         EVXLIQ(N,NP)=0.d0
2 	   EVX(N,NP)=0.D0	
      
       DO 3 NG=1,NPOL
	  X=R2(NG)
	  W=W2(NG)
	  E=TKI*X

c******************************************************************************
c!!!!	POINTS TO NOTE: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c	
c	 1. To introduce new elastic cross-sections you need only modify 
c		between the two lines of ************
c	 2. For real gases comment out all references to sy0 and syl
c		in the model section and the converse for model gases
c	 3.	Always set both sy0 and syl:
c			- For isotropic scattering set syl=0.d0
c			- For anisotropic scattering syl must be specified
c	  					
c------------------------------------------------------------------------------
									! REAL GASES
	  SY0=SIGMA(JJJ,0,E)
c	  SYL=SIGMA(JJJ,L,E)	! For anisotropic scattering
		
	  SYL=0.d0				! For isotropic scattering

c------------------------------------------------------------------------------
									! MODEL GASES

c..............................................................................
                                    ! Hg-like gas
c        qval = 1.6d0
c        sigma0val = 280.d0
c        x1val = 10.d0*(e-0.5d0)
c        x2val = 0.5d0*(e-0.8d0)
             
c        a02angstromsq=((5.291772108d-11)**2.d0)/1d-20        ! 1 a0 to m
c        sy0=(qval+x1val)*(qval+x1val)/(1.d0+x1val*x1val)
c     $  /(1.d0+x2val*x2val)*sigma0val+10.d0 
c        sy0=sy0*a02angstromsq
c        syl=0.d0       
c        write(99,*) e,sy0
c..............................................................................

c	  sy0=6.d0/dsqrt(e)           ! Maxwell model
c	  if (e.eq.0.d0) sy0=100.d0
c	  syl=0.d0
c..............................................................................

c	  sy0=0.d-4					  ! Hard sphere model
c	  syl=0.d0

c	  sy0=10.d0
c	  syl=0.d0
c..............................................................................

c	  sy0=10.d0							! Reid Model B
c	  syl=0.d0
c	  if (l.eq.2) syl=5.714285714d0
c	  if (l.eq.4) syl=1.26984127d0

c..............................................................................

										! Reid Model C	     
c		tmp=10.97676d0
c		sy0=tmp*0.63347d0
c		if (l.eq.0) syl=sy0
c		if (l.eq.1) syl=tmp*(-0.277541d0)
c		if (l.eq.2) syl=tmp*(0.0783929d0)
c		if (l.eq.3) syl=tmp*(-0.0162315d0)
c		if (l.eq.4) syl=tmp*(0.00264603d0)
c		if (l.eq.5) syl=tmp*(-0.000355295d0)
c		if (l.eq.6) syl=tmp*(0.000040532d0)
c		if (l.eq.7) syl=tmp*(-4.01799d-6)
c		if (l.eq.8) syl=tmp*(3.52087d-7)	
c		if (l.eq.9) syl=tmp*(-2,76486d-8)
c		if (l.eq.10) syl=tmp*(1.62981d-9)
		   
c..............................................................................
		
c		   tmp=18.29141378d0			! Reid Model D
c		   sy0=tmp*0.3802016314099497d0
c		   if (l.eq.0) syl=sy0
c		   if (l.eq.1) syl=tmp*(-0.16650296685856297)
c		   if (l.eq.2) syl=tmp*0.252983735017978
c		   if (l.eq.3) syl=tmp*(-0.027682353217118552)
c		   if (l.eq.4) syl=tmp*0.0760334106091381
c		   if (l.eq.5) syl=tmp*0.06424451994922523
c		   if (l.eq.6) syl=tmp*(-0.010964021071413504)
c		   if (l.eq.7) syl=tmp*0.045372021147647094
c		   if (l.eq.8) syl=tmp*(-0.004082795114054207)
c		   if (l.eq.8) syl=tmp*(-0.01242501356770807)
c	 	   if (l.eq.10) syl=tmp*0.010362738908258541
c
c..............................................................................
					
									! Double ramp model
c	      sy0=10.d0
c		  syl=0.d0
		
c------------------------------------------------------------------------------
! GASES:

!	  SY0=SIGMA(JJJ,0,E)
!	  SYL=SIGMA(JJJ,L,E)	! For anisotropic scattering
!!        SYL=0.d0            ! Isotropic scattering

! LIQUIDS:  
!        sy0bar=SIGMA(JJJ,0,E)   !
!        sylbar=SIGMA(JJJ,L,E)    ! For isotropic with structure

!        sy0=SIGMA(JJJ+1,0,E)    ! Gas form of elastic cross-section
!c        write(509,*) e,jjj,sy0bar,sylbar,sy0
!        syl=0.d0    ! For isotropic
c******************************************************************************

!        sybar=sy0bar-sylbar
        sy=sy0-syl  
	  if (l.eq.0) then 
	      sybar=0.d0
	      sy=0.d0
	  end if
		  
       
          CALL SONP(AL,NLIM,X,S)
          DO 4 N=1,NLIM
	     if (LL.eq.1) then
		    ! Interaction integral for process jjj for rate calculations
	       vrate0(no_gas,0,n-1)=vrate0(no_gas,0,n-1)+w*x**(fl+1.d0)*
     $			s(n)*s(1)*sy0*sca2
		   vrate1(no_gas,0,n-1)=vrate1(no_gas,0,n-1)+w*x**(fl+1.d0)*
     $			s(n)*s(2)*sy0*sca2
	       	! Calculating the direct part
	     endif
         DO 4 NP=1,N
	  if (LL.GT.1) THEN
            EVX(N,NP)=EVX(N,NP)+W*X**(FL+1.D0)*S(N)*S(NP)*SY
            EVXLIQ(N,NP)=EVXLIQ(N,NP)+W*X**(FL+1.D0)*S(N)*S(NP)*SYBAR
	    EVX(NP,N)=EVX(N,NP)    
	    EVXLIQ(NP,N)=EVXLIQ(N,NP)    
			! Interaction integrals are symmetric in nu indicies
	  endif
4 	  CONTINUE
3      CONTINUE
6     DO 512 N=1,NLIM
        DO 5 NP=1,NLIM
         VXELAS(LL,N,NP)=EVX(N,NP)*SCA2
         VXLIQ(LL,N,NP)=EVXLIQ(N,NP)*SCA2
5       CONTINUE
512	CONTINUE

11    CONTINUE

c	Output formatting

505	FORMAT(1X,'L = ',I2)
	WRITE(6,998)
998	FORMAT(' I AM AT VINT END')
	WRITE(2,2505)NPOL,WI
2505	FORMAT(' FOR ELASTICS TRAPEZOIDAL RULE USED ; NO. OF PTS = '
     1  ,I4,' WIDTH = ',F8.6)

      RETURN
	END


c ****************************************************************************

       SUBROUTINE UINT(LLIM,NLIM,TKI,TK,NGQ1,NPOL,WI,ix1,ix2,NO_GAS)

c	This subroutine calculates the sw*v's for the inelastic processes
c 	starting at the cross-section ix1 and finishing at ix2
c	i.e. exp(-E_i/kTb)/Z_o*v's

c	Ensure inelastic cross-sections are between ix1 and ix2

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/WEALTH/FACLOG(150),GAMLOG(150),VX(10,100,100),A(100,100,9),
     +   VXLIQ(10,100,100),VXATT(10,100,100),VXELAS(10,100,100),
     +      VXINEL(10,100,100)      
      COMMON/LAW/SW(21),DE(21),V(21),VI(21)
      COMMON/GOOD/EN(21,1100),QS(21,1100,0:10),NO(21),npq(21)
      COMMON/INLAW/nq,nq1,nq2,nqi1,nqi2
	COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)
      COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)
	DIMENSION SKJ(100)
      
	if(ix2.lt.0)then
	 print*,'No inelastics for this gas'
	return
	end if

	print*, 'In Uint'

	ALF=0.0D0
	SCA1=1.D0
	SCA2=1.D-26		! Scaling from cotes
	CALL COTES(NPOL,WI,1)
	

	do 103 j=1,30
	 do 103 n=1,nlim
	   vrate0(no_gas,j,n-1)=0.d0
	   vrate1(no_gas,j,n-1)=0.d0
	   vrate0(no_gas,-j,n-1)=0.d0
	   vrate1(no_gas,-j,n-1)=0.d0
103	 continue


	DO 11 LL=1,LLIM		! l-loop
          L=LL-1
	  WRITE(6,521)L
          FL=DFLOAT(L)
          AL=FL+0.5D0
          FL2=0.5D0*FL
     
          DO 2 N=1,NLIM
           DO 2 NP=1,NLIM
2           EVX(N,NP)=0.D0
     
        DO 3 NG=1,NPOL		! Energy integration
	   X=R2(NG)
	   W=W2(NG)
         CALL SONP(AL,NLIM,X,S)
         IF(X.EQ.0.D0)THEN
          DUM1=1.D0
          DUM2=1.D0
         ELSE
          DUM1=X**FL2
          DUM2=X**ALF
         ENDIF


      no_process=0		! Specifies the process number for rate calculations
						! Elastics: no_process=0; 
						! Inelastics etc: no_process>0;
						! Superelastics: no_process=-no_process for inelastics.

	   
 	   DO 4 JA=ix1,ix2	! Process loop
	      
		  no_process=no_process+1

		  JAY=JA
            XKJ=X+DE(JAY)/TKI
            E=XKJ*TKI
            XKJFL2=XKJ**FL2
            CALL SONP(AL,NLIM,XKJ,SKJ)

c******************************************************************************
c!!!!	POINTS TO NOTE: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c	
c	 1. To introduce new inelastic cross-sections you need only modify 
c		between the two lines of ************
c	 2. For real gases comment out all references to sig0 and sigl
c		in the model section and the converse for model gases
c	 3.	Always set both sig0 and sigl:
c			- For isotropic scattering set sigl=0.d0
c			- For anisotropic scattering sigl must be specified

c------------------------------------------------------------------------------
										! REAL GASES
c	    SIG0=SIGMA(JAY,0,E)
c	    SIGL=SIGMA(JAY,L,E)		! For anisotropic scattering

c 		SIGL=0.d0				! For isotropic scattering

c------------------------------------------------------------------------------
										! MODEL GASES
c..............................................................................

c	    if (e.ge.0.2d0) then        ! Reid ramp model
c	     sig0=10.0d0*(e-0.2d0)
c	    else
c	     sig0=0.d0
c	    endif
c		sigl=0.d0
c..............................................................................
				
c		if (e.ge.0.5d0) then 
c		  sig0=0.1d0
c		else
c		  sig0=0.d0
c		endif
c		sigl=0.d0
		
c..............................................................................

c	    if (e.ge.0.516d0) then		! H2 model 
c	     sig0=4.0d-1*(e-0.516d0)
c	    else
c	     sig0=0.d0
c	    endif
c		sigl=0.d0

c..............................................................................
		
c		if (ja.eq.ix1)  then		! H2 swarm/theory	
c	      sig0=0.d0
c	      if (e.ge.0.044d0) then
c			sig0=0.5d0
c	      endif
c		elseif (ja.eq.ix2) then
c	      sig0=0.d0
c		  if (e.ge.0.516d0) then
	        
c			slope=0.15d0	! Swarm
c			slope=0.25d0	! Theory	
c							! Ang^2 per eV
c		    enstar=0.516d0+0.5d0/slope
c		    
c			if (e.le.enstar) then
c		     sig0=slope*(e-0.516d0)
c			else 
c			 sig0=0.5d0
c			endif
!c		  write(56,*) ja,e,sig0
c	      endif
c		 endif
c	     sigl=0.d0  ! Isotropic scattering   

c..............................................................................
	     
c		 sig0=0.d0					! Reid anisotropic B model		 
c		 sigl=0.d0
c		 if (e.ge.0.516d0) then
c	       sig0=0.4d0*(e-0.516d0)
c	       if (l.eq.0) sigl=sig0
c	       if (l.eq.2) sigl=0.2285714286d0*(e-0.516d0)
c		   if (l.eq.4) sigl=0.05079365079d0*(e-0.516d0)
c		 endif

c..............................................................................
		
		 							! Reid Model C	     
c		 sig0=0.d0					
c		 sigl=0.d0
c		 tmp=0.631443d0*(e-0.516d0)
c		 sig0=tmp*0.63347d0
c		 if (e.ge.0.516d0) then 
c
c		if (l.eq.0) sigl=sig0
c		if (l.eq.1) sigl=tmp*(-0.277541d0)
c		if (l.eq.2) sigl=tmp*(0.0783929d0)
c		if (l.eq.3) sigl=tmp*(-0.0162315d0)
c		if (l.eq.4) sigl=tmp*(0.00264603d0)
c		if (l.eq.5) sigl=tmp*(-0.000355295d0)
c		if (l.eq.6) sigl=tmp*(0.000040532d0)
c		if (l.eq.7) sigl=tmp*(-4.01799d-6)
c		if (l.eq.8) sigl=tmp*(3.52087d-7)	
c		if (l.eq.9) sigl=tmp*(-2,76486d-8)
c		if (l.eq.10) sigl=tmp*(1.62981d-9)
c
c	     endif

		 
c..............................................................................
		 
c		 sig0=0.d0					! Reid anisotropic D model		 
c		 sigl=0.d0
c		 tmp=1.052073343d0*(e-0.516d0)
c		 if (e.ge.0.516d0) then
c	       sig0=tmp*0.3802016314099497d0
c		   if (l.eq.0) sigl=sig0
c		   if (l.eq.1) sigl=tmp*(-0.16650296685856297)
c		   if (l.eq.2) sigl=tmp*0.252983735017978
c		   if (l.eq.3) sigl=tmp*(-0.027682353217118552)
c		   if (l.eq.4) sigl=tmp*0.0760334106091381
c		   if (l.eq.5) sigl=tmp*0.06424451994922523
c		   if (l.eq.6) sigl=tmp*(-0.010964021071413504)
c		   if (l.eq.7) sigl=tmp*0.045372021147647094
c		   if (l.eq.8) sigl=tmp*(-0.004082795114054207)
c		   if (l.eq.8) sigl=tmp*(-0.01242501356770807)
c	 	   if (l.eq.10) sigl=tmp*0.010362738908258541
c		 endif
		 
c..............................................................................

									! Hg model
				
		if (e.ge.(4.9d0)) then 
		  sig0=2.d0
		else
		  sig0=0.d0
		endif
		sigl=0.d0	

c------------------------------------------------------------------------------
c******************************************************************************
		
		  if (l.eq.0) sigl=sig0  
            write(734,*) e,sig0

 	      B1=VI(JAY)*XKJ**(FL2+1.D0)
            B2=SIG0*XKJFL2
            B3=SIGL*DUM1
            C1=V(JAY)*DUM1*XKJ
            C2=SIG0*DUM1
            C3=SIGL*XKJFL2
            C4=SW(JAY)

	     DO 41 N=1,NLIM	! nu-loops
	       if (l.eq.0) then ! Only for l.eq.0
	
      vrate0(no_gas,no_process,n-1)=vrate0(no_gas,no_process,n-1)+
     +			(C4/DUM2)*B1*SKJ(N)*W*B2*SKJ(1)*sca2
	vrate0(no_gas,-no_process,n-1)=vrate0(no_gas,-no_process,n-1)+
     +						(C4/DUM2)*C1*S(1)*W*(C2*S(N))*sca2
	vrate1(no_gas,no_process,n-1)=vrate1(no_gas,no_process,n-1)+
     +						(C4/DUM2)*B1*SKJ(N)*W*B2*SKJ(2)*sca2
	vrate1(no_gas,-no_process,n-1)=vrate1(no_gas,-no_process,n-1)+
     +						(C4/DUM2)*C1*S(2)*W*(C2*S(N))*sca2
			! For rate calculations exp(-E_i/kTb)*v(0,0,nud)direct
			! NB It has to be like this to get each process	
			! individually.  Note the integration loop is outside the process loop. 
	       endif
              DO 410 NP=1,NLIM
                  EVX(N,NP)=EVX(N,NP)+(C4/DUM2)*(B1*SKJ(NP)*W*(B2*SKJ(N)
     1                -B3*S(N))+C1*S(NP)*W*(C2*S(N)-C3*SKJ(N)))

410          CONTINUE
41	    CONTINUE
	
4	 CONTINUE   ! End of process JA loop  
3     CONTINUE	! End of l loop

      DO 6 N=1,NLIM
         DO 6 NP=1,NLIM
	        VXINEL(LL,N,NP)=SCA2*EVX(N,NP)+VXINEL(LL,N,NP)
6	CONTINUE

11    CONTINUE	! End of l-loop

c	Outputs

521	FORMAT(1X,'L = ',I3)
321	FORMAT(' START GAUSS. QUAD. ')
322	FORMAT(1X,I5,F8.4)
432	FORMAT(' END GAUSS. QUAD. ')
543	FORMAT(' START INT. ')
654	FORMAT(' END INT. ')
	WRITE(2,2505)NPOL,WI
2505	FORMAT(' FOR INELASTICS SIMPSONS RULE USED ; NUMBER OF
     1 POINTS = ',I4,' WIDTH = ',F8.6)

      RETURN
	END

c ****************************************************************************

      SUBROUTINE ATTINT(LLIM,NLIM,TKI,NGQ1,NPOL,WI,JJJ,NO_GAS)

c	This subroutine calculates the v's for attachment process
c    	which is defined by the cross-section jjj.

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/WEALTH/FACLOG(150),GAMLOG(150),VX(10,100,100),A(100,100,9),
     +   VXLIQ(10,100,100),VXATT(10,100,100),VXELAS(10,100,100),
     +      VXINEL(10,100,100)
      COMMON/GOOD/EN(21,1100),QS(21,1100,0:10),NO(21),npq(21)
      COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)
    
      COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)
      
      
      WRITE(6,123)NPOL,WI
 123	FORMAT(' I AM AT ATTINT ',I6,F8.5)

	CALL COTES(NPOL,WI,1)
	AA=0.D0
	SCA1=1.D0
	SCA2=1.D-26
	NGQ=NPOL
	
	DO 11 LL=1,LLIM   ! l-loop
          L=LL-1
	  WRITE(6,505)L
505	  FORMAT(1X,'L = ',I2)
          FL=DFLOAT(L)
	  AL=FL+.5D0
          DO 2 N=1,NLIM
           DO 2 NP=1,NLIM
2           EVX(N,NP)=0.D0
      
	 DO 3 NG=1,NGQ  	! Energy integration
	  X=R2(NG)
	  W=W2(NG)
	  E=TKI*X
 	  SY= SIGMA(JJJ,0,E)
        CALL SONP(AL,NLIM,X,S)
        DO 4 N=1,NLIM	! nu-loops
         DO 4 NP=1,N
           EVX(N,NP)=EVX(N,NP)+W*X**(FL+1.D0-AA)*S(N)*S(NP)*SCA1*SY
	     EVX(NP,N)=EVX(N,NP)
4	  CONTINUE
3      CONTINUE
	
	print*, 'nqatt= ',jjj
6	 DO 512 N=1,NLIM
          DO 5 NP=1,NLIM
         VXATT(LL,N,NP)=VXATT(LL,N,NP)+EVX(N,NP)*SCA2 
 5 	 CONTINUE  !VX(LL,N,NP)=VX(LL,N,NP)+EVX(N,NP)*SCA2
 	    vrate0(no_gas,jjj,n-1)=vx(1,1,n)   ! Rate coefficient
				!**** BE CAREFUL FOR GAS 2 in CALCULATING RATES
512    continue

11     CONTINUE		! End of l-loop

	 WRITE(6,998)
998	 FORMAT(' I AM AT ATTINT END')
	 WRITE(2,2505)NPOL,WI
2505	 FORMAT(' FOR attachment TRAPEZOIDAL RULE USED ; NO. OF PTS = '
     1  ,I4,' WIDTH = ',F8.6)
       
	 RETURN
	 END

c ****************************************************************************

      SUBROUTINE RINT(LLIM,NLIM,TKI,NGQ1,NPOL,WIDTH,DELTA,IOP,
     +											NQION,NO_GAS)

c	This subroutine calculates the v's for the ionization process NQION
c	and stores them in vxion (which is to zeroth order in the mass ratio)

      IMPLICIT DOUBLE PRECISION(a-h,o-z) 
      COMMON/ION/AA,SIG(5000)
      COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)
      COMMON/COL/VXION(10,100,100)
      COMMON/LAW/SW(21),DE(21),V(21),VI(21)
      COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)

	WRITE(6,999)NPOL,WIDTH
999	FORMAT(' I AM AT RINT',I6,F8.5)

	XI=DE(NQION)/TKI
	SW(NQION)=1.d0
	NGQ=NPOL
	AA=0.0d0
	SCA1=1.d0
	SCA2=1.d-26
	CALL COTES(NPOL,WIDTH,2)
43	DO 10 I=1,NGQ
 	 E=(R2(I)+XI)*TKI
10	 SIG(I)=SIGMA(NQION,0,E)
       
	 DO 11 LL=1,LLIM		! l-loop

          L=LL-1
 	  WRITE(6,521)L
521	  FORMAT(1X,'L = ',I3)
          FL=L
          AL=FL+0.5d0
	  FL2=FL+1.d0
       
	  DO 2 N=1,NLIM
           DO 2 NP=1,NLIM
2           EVX(N,NP)=0.d0
       
	    DO 3 NG=1,NGQ		! Energy integration

	     X=R2(NG)
	     W=W2(NG)
   	     XXI=X+XI
             CALL SONP(AL,NLIM,XXI,S)
   	     SIGT=SIG(NG)
	     B1=VI(NQION)*XXI**FL2*SIGT
	     if(X.eq.0.d0)then
	      XAA = 1.d0
	     else
	      XAA = X**AA
	     end if

        DO 41 N=1,NLIM		! nu-loop
         DO 41 NP=1,N
	    EVX(N,NP)=EVX(N,NP)+SW(NQION)/XAA*B1*S(NP)*W*S(N)*SCA1
41      CONTINUE

3      CONTINUE			! End of Energy integration

	 DO 5 N=1,NLIM
 	  DO 5 NP=1,N
5	   EVX(NP,N)=EVX(N,NP)  ! Symmetric direct part
       
	 DO 6 N=1,NLIM
          DO 6 NP=1,NLIM
6          VXION(LL,N,NP)=SCA2*EVX(N,NP)
         
11     CONTINUE			! l-loop

	print*, 'nqion= ',nqion
	 DO 61 N=1,NLIM
	    vrate0(no_gas,nqion,n-1)=vxion(1,1,n)    ! For rate calculations
					   ! This is the direct part of the integral 	
61	 CONTINUE

	print*,' At Call REST'
	CALL REST(NLIM,NGQ,TKI,XI,DELTA,IOP,NQION)
	print*,' End Call REST'

C	IF IOP = 1 , THEN ALL FRACTIONS ARE EQUI PROBABLE

	DO 7 N=1,NLIM
	 DO 7 NP=1,NLIM
7	 VXION(1,N,NP)=VXION(1,N,NP)-2.d0*SW(NQION)*EVX(N,NP)



c	Output formatting

	WRITE(2,2505)NPOL,WIDTH,DELTA,IOP
2505	FORMAT(' FOR IONIZATION SIMPSONS RULE USED ; NUMBER OF POINTS
     1 = ',I4,' WIDTH = ',F8.4/' DELTA = ',F8.4,' IOP = ',I2)

2909	RETURN
	END

	SUBROUTINE REST(NLIM,NGQ,TKI,XI,DELTA,IOP,NQION)

C	CALCULATES RESTITUTING PART OF L = 0 INTERACTION INTEGRAL

	IMPLICIT DOUBLE PRECISION(a-h,o-z) 
	DIMENSION SK(5000),SJKI(5000),SD(5000)
	DIMENSION GK(5000)
	COMMON/ION/AA,SIG(5000)
	COMMON/LAW/SW(21),DE(21),V(21),VI(21)
	COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)

	SCA1=1.d-18
	SCA2=1.d-8

	DO 5 N=1,NLIM
 	 GK(N)=0.d0
	 DO 5 NP=1,NLIM
 5	  EVX(N,NP)=0.d0
	
	IF(IOP.NE.1)GO TO 50
	
	DO 30 K=1,NGQ  	 ! First integration
	 
	 XK=R2(K)
	 WK=W2(K)
	 if(XK.eq.0.d0)then
	  XKAA = 1.d0
	 else
	  XKAA = XK**AA
	 end if
	 CALL SONP(.5d0,NLIM,XK,SK)
	 
	 DO 15 J=1,NGQ  	! Second integration
	
	  XJ=R2(J)
	  WJ=W2(J)
	  XJK=XJ+XK
	  XJKI=XJK+XI
	  CALL SONP(.5d0,NLIM,XJKI,SJKI)
 	  E=XJKI*TKI
	  SIGT=SIGMA(NQION,0,E)
	  if(XJ.eq.0.d0)then
	   XJAA = 1.d0
	  else
	   XJAA = XJ**AA
	  end if
	  if(XJK.eq.0.d0)then
	   XJKIXJK = 0.d0
	  else
	   XJKIXJK = XJKI/XJK
	  end if
	  B1=XJKIXJK*SIGT*SCA1/XJAA
	
	  DO 10 NP=1,NLIM
 10	   GK(NP)=GK(NP)+B1*SJKI(NP)*WJ

 15	 CONTINUE    	! End of second energy integration
 
	 DO 20 NP=1,NLIM
 20	  GK(NP)=GK(NP)*SCA2/XKAA
	
	 DO 25 N=1,NLIM
	  DO 25 NP=1,NLIM
 25	   EVX(N,NP)=EVX(N,NP)+WK*SK(N)*GK(NP)*SCA1

 30	CONTINUE		!  End of first integration

	DO 35 N=1,NLIM
	 DO 35 NP=1,NLIM
 35	  EVX(N,NP)=EVX(N,NP)*SCA2*VI(NQION)

	RETURN  		! If IOP.eq.1

 50	CONTINUE		! If IOP.ne.1

	DELTA2=1-DELTA

	DO 70 K=1,NGQ	! Energy integration
	 W=W2(K)
	 X=R2(K)
	 XXI=X+XI
	 if(X.eq.0.d0)then
	  XAA = 1.d0
	 else
	  XAA = X**AA
	 end if
	 DEL1X=X*DELTA
	 DEL2X=X*DELTA2
	 CALL SONP(.5d0,NLIM,XXI,SK)
	 CALL SONP(.5d0,NLIM,DEL1X,SJKI)
	 CALL SONP(.5d0,NLIM,DEL2X,SD)
	 SIGT=SIG(K)
	 B1=XXI*SIGT*SCA1/XAA
	 DO 60 N=1,NLIM
	  DO 60 NP=1,NLIM
 60	   EVX(N,NP)=EVX(N,NP)+B1*SK(NP)*W*(SJKI(N)+SD(N))
 70	CONTINUE
	
	DO 65 N=1,NLIM
	 DO 65 NP=1,NLIM
 65	  EVX(N,NP)=SCA2*EVX(N,NP)*VI(NQION)/2.d0

	RETURN
	END

C	**********************************************************

      SUBROUTINE COLMAT(WM,C,LMAX,RMAX,T,TI,NO_GAS)

c	This routine calculates the dimensionless V's from the 
c        dimensionless v's (and sw*v's for 
c	inelastic processes) and stores them in VX and vxion.
c	These are then converted into the dimensionless J's stored
c	in A.  
c	These J's are then converted into dimensionless a's
c	in the calling routine.



      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER R,S,RMAX
      COMMON/WEALTH/FACLOG(150),GAMLOG(150),VX(10,100,100),A(100,100,9),
     +   VXLIQ(10,100,100),VXATT(10,100,100),VXELAS(10,100,100),
     +      VXINEL(10,100,100)
      common/col/vxion(10,100,100)
      COMMON/RATES/vrate0(2,-30:30,0:100),vrate1(2,-30:30,0:100)
      COMMON/INLAW/nq,nq1,nq2,nqi1,nqi2

	WRITE(6,999)
999	FORMAT(' I AM AT COLMAT	')

	NLIM=RMAX+1		! For first order massrat only
	LLIM=LMAX+1		! For first order massrat only
	
c	Converting the v's and sw*v's to V's and sw*V's

	DO 5 LL=1,LLIM
	 DO 5 N=1,NLIM
	  AF=FACLOG(N)-GAMLOG(LL+N)
	  DO 5 NP=1,NLIM
	   AFP=FACLOG(NP)-GAMLOG(LL+NP)
	   AL=DEXP(AF+AFP)
	   SQAL=DSQRT(AL)
	   vxion(ll,n,np)=sqal*vxion(ll,n,np)
	   VXELAS(ll,n,np)=sqal*VXELAS(ll,n,np)
	   VXINEL(ll,n,np)=sqal*VXINEL(ll,n,np)
	   VXATT(ll,n,np)=sqal*VXATT(ll,n,np)
	   VXLIQ(ll,n,np)=sqal*VXLIQ(ll,n,np)
5	   VX(LL,N,NP)=SQAL*VX(LL,N,NP)

	DO 6 I=-NQ,NQ
	 AF=FACLOG(1)-GAMLOG(2)
	 DO 6 NP=1,NLIM
	   AFP=FACLOG(NP)-GAMLOG(1+NP)
	   AL=DEXP(AF+AFP)
	   SQAL=DSQRT(AL)
	    DO 6 NO_GAS=1,2
		vrate0(no_gas,i,np-1)=vrate0(no_gas,i,np-1)*sqal
6	CONTINUE

	DO 7 I=-NQ,NQ
	 AF=FACLOG(2)-GAMLOG(3)
	 DO 7 NP=1,NLIM
	   AFP=FACLOG(NP)-GAMLOG(1+NP)
	   AL=DEXP(AF+AFP)
	   SQAL=DSQRT(AL)
	    DO 7 NO_GAS=1,2
	    vrate1(no_gas,i,np-1)=vrate1(no_gas,i,np-1)*sqal
7	CONTINUE

c	Mass ratio expansion to first order in mass ratio.
c	These are the dimensionless J's

      DO 1 L=1,LMAX
       FL=DFLOAT(L-1)
       DO 1 R=1,RMAX
        FR=DFLOAT(R-1)
        TR=FR+FL+.5D0
        DO 1 S=1,RMAX
         FS=DFLOAT(S-1)
         TS=FS+FL+.5D0
!         TERM1=VXLIQ(L,R,S)+VXATT(L,R,S)+VXINEL(L,R,S)+vxion(l,r,s)     ! Liquids
         TERM1=VXELAS(L,R,S)+VXATT(L,R,S)+VXINEL(L,R,S)+vxion(l,r,s)  ! Gases
         
         TERM2=0.D0
         TERM3=0.d0
         
         if (l.eq.1) then  ! THE L=0 operator
         IF(R.GT.1)TERM2=TERM2+DSQRT(FR*TR)*VXELAS(L,R-1,S)
         TERM2=TERM2+DSQRT((FS+1.D0)*(TS+1.D0))*VXELAS(L,R,S+1)
         IF(R.GT.1)TERM2=TERM2-2.D0*(FL+1.D0)/(2.D0*FL+1.D0)
     1    *DSQRT(FR*(TS+1.D0))*VXELAS(L+1,R-1,S)
         IF(L.GT.1)TERM2=TERM2-2.D0*FL/(2.D0*FL+1)*
     1    DSQRT((FS+1.)*TR)*VXELAS(L-1,R,S+1)
         TERM3=0.D0
         IF(R.GT.1.AND.S.GT.1)TERM3=TERM3+(FL+1.D0)/(2.D0*FL+1.D0)
     1    *2.D0*DSQRT(FR*FS)*VXELAS(L+1,R-1,S-1)
         IF(L.GT.1)TERM3=TERM3+FL/(2.D0*FL+1.D0)*2.D0*DSQRT(TR*TS)
     1    *VXELAS(L-1,R,S)
         TERM3=TERM3-(FL+FR+FS)*VXELAS(L,R,S)
         end if
         
         TERM=(TERM1+WM*((1.D0-T/TI)*TERM2+TERM3))
         A(R,S,L)=A(R,S,L)+TERM*C
1     CONTINUE
	WRITE(6,99)
99	FORMAT(' AT COLMAT END')
      RETURN
	END

c ****************************************************************************

       FUNCTION SIGMA(JJ,L,E)

       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       COMMON/GOOD/EN(21,1100),QS(21,1100,0:10),NO(21),npq(21)
       COMMON/LAW/SW(21),DE(21),V(21),VI(21)
       COMMON/INLAW/nq,nq1,nq2,nqi1,nqi2 
       DIMENSION X(1100),Y(1100)
C      WRITE(3,2001)NQ,JJ
 2001  FORMAT(1X,I4,5X,I4)
	j=jj

!	if(j.eq.nq.and.e.lt.1.0d0)then   !NB ONLY VALID FOR AR
!	  z=phase(e)
!	else
        IMAX=NO(J)
	 DO 1 I=1,IMAX
	 X(I)=EN(J,I)
	 Y(I)=QS(J,I,L)
   1	 CONTINUE
	 sigma=FLAGR(X,Y,E,1,IMAX)
	 
!	end if
 10	RETURN
	END

	function phase(en)
	implicit double precision(a-h,o-z)
C
C  PARAMETERS OF PHASE SHIFT ANALYSIS. From Steve Biagi's 1994 data file.
C  ###### For ARGON ########
      PIR2=8.79735669D-17
      ARY=13.6056981
      APOL=11.08 
      LMAX=100 
      AA=-1.488 
      DD=65.4
      FF=-84.3 
      E1=0.883 
      API=3.141592654 
	phase = 0.0d0
C 
      IF(EN.GT.1.0) GO TO 100
      IF(EN.EQ.0.0) phase=7.79E-16*1.0d16
      IF(EN.EQ.0.0) GO TO 100 
      AK=DSQRT(EN/ARY) 
      AN0=-AA*AK*(1.0+(4.0*APOL/3.0)*AK*AK*DLOG(AK))-(API*APOL/3.0)*AK*A
     /K+DD*AK*AK*AK+FF*AK*AK*AK*AK 
      AN1=(API/15.0)*APOL*AK*AK*(1.0-DSQRT(EN/E1))
      AN2=API*APOL*AK*AK/105.0 
      AN0=DATAN(AN0) 
      AN1=DATAN(AN1)
      AN2=DATAN(AN2)
      SUM=(DSIN(AN0-AN1))**2
      SUM=SUM+2.0*(DSIN(AN1-AN2))**2
      DO 10 J=2,LMAX-1 
      SUMI=6.0/((2.0*J+5.0)*(2.0*J+3.0)*(2.0*J+1.0)*(2.0*J-1.0))
      SUM=SUM+(J+1.0)*(DSIN(DATAN(API*APOL*AK*AK*SUMI)))**2 
   10 CONTINUE 
      phase=SUM*4.0*PIR2/(AK*AK)*1.0d16
 100	continue
	return
	end

c ****************************************************************************

      FUNCTION FLAGR(X,Y,XARG,IDEG,N)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(1100),Y(1100)

C     FINDS THE POSITION OF XARG WITH RESPECT TO THE ARRAY X(I)

      IF(XARG.LE.X(2))GO TO 100
      IF(XARG.GE.X(N-1))GO TO 101
      NA=0
      NB=N
  10  K=(NA+NB)/2
      IF(XARG-X(K))1,2,3
   1  NB=K
      IF(NB-NA-1)4,4,10
   3  NA=K
      IF(NB-NA-1)4,4,10
   2  FLAGR=Y(K)
      RETURN
   4  MIN=NA-IDEG/2
C
C     LAGRANGIAN INTERPOLATION
C
      FACTOR=1.D0
      MAX=MIN+IDEG
      J=MIN
   20 IF(XARG.NE.X(J))GO TO 21
      FLAGR=Y(J)
      RETURN
  21  FACTOR=FACTOR*(XARG-X(J))
      J=J+1
      IF(J.LE.MAX)GO TO 20
      YEST=0.D0
      DO 25 I=MIN,MAX
      TERM=Y(I)*FACTOR/(XARG-X(I))
      DO 24 J=MIN,MAX
      IF(I.NE.J)TERM=TERM/(X(I)-X(J))
  24  CONTINUE
  25  YEST=YEST+TERM
      FLAGR=YEST
      RETURN
 100  FLAGR=Y(1)+(Y(2)-Y(1))*(XARG-X(1))/(X(2)-X(1))
      RETURN
 101  FLAGR=Y(N)+(Y(N-1)-Y(N))*(XARG-X(N))/(X(N-1)-X(N))
      RETURN
	END

c ****************************************************************************

      SUBROUTINE SONP(AL,NLIM,X,S)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c      DIMENSION S(62)
       DIMENSION S(100)
      S(1)=1.D0
      IF(NLIM.EQ.1)GO TO 10
      S(2)=AL+1.D0-X
      IF(NLIM.EQ.2)GO TO 10
      DO 1 NNU=3,NLIM
      FK=NNU-1
   1  S(NNU)=((2.D0*FK+AL-1.D0-X)*S(NNU-1)-(FK+AL-1.D0)*S(NNU-2))/FK
   10 RETURN
	END

c ****************************************************************************

	SUBROUTINE COTES(NGQ,WIDTH,IFACT)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	COMMON S(100),R2(5000),W2(5000),EVX(100,100),EVXLIQ(100,100),
     + EVXATT(100,100)
	V1=4.D0/3.D0
	V2=2.D0/3.D0
  	H=WIDTH/DFLOAT(NGQ-1)
	ALOGH=DLOG(H)
      SCLOG=DLOG(1.D26)
  	DO 100 NG=1,NGQ
 	Y=DFLOAT(NG-1)*H
	X=Y/(1.D0-Y)
	W=-X-2.D0*DLOG(1.D0-Y)+ALOGH+SCLOG
 	W=DEXP(W)
	R2(NG)=X
	GO TO (10,15),IFACT
C ***** TRAPEZOIDAL RULE
 10	W2(NG)=W
   	IF(NG.EQ.1.OR.NG.EQ.NGQ)W2(NG)=W/2.D0
	GO TO 100
C ***** SIMPSON'S RULE
 15	KK=NG/2
	KK=2*KK
	W2(NG)=W*V1
	IF(KK.NE.NG)W2(NG)=W*V2
	IF(NG.EQ.1.OR.NG.EQ.NGQ)W2(NG)=W/3.D0
 100	CONTINUE
	RETURN
	END

